% variance decomposition of detrended (or demeaned) output, inflation, and
% interest rate 


% Taeyoung Doh
% created: 12/17/2010
%--------------------------------------------------------------------1.0
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;


%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);




global datasel psel 
% MHRUN
%  length of mhrun defines the number of columns in graph table
%	odd integer  (e.g., 1) for linear posterior draws
%	even integer (e.g., 2) for quadratic posterior draws
 
mhrun = 1;
mhrunid  = num2str(mhrun','%2.0f');

block1   = 11;		% starting block
blockn   = 150;		% ending block
nsim     = 10000;	% number of simulation draws in each block

hpdprob  = 0.90;	% probability for highest posterior density


%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% parameter names
pnames = feval(fnames_.pnames,msel);

% reset diary file
diaryfile = [path_.result 'stat_mhrun' mhrunid '.log'];
if exist(diaryfile,'file')
	delete(diaryfile);
end
diary(diaryfile);
diary off;

%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------


% load path
  lpath = [path_.para 'mhrunpm' mhrunid '/' ];  
% initialize output
nsimul = nsim/100;
        Vout = zeros(4,3,nsimul*(blockn-block1+1));   drawout = zeros(2,4,3); doutmean=zeros(4,3);
        Vinf = zeros(4,3,nsimul*(blockn-block1+1));   drawinf = zeros(2,4,3); dinfmean=zeros(4,3);
        Vint = zeros(4,3,nsimul*(blockn-block1+1));   drawint = zeros(2,4,3); dintmean=zeros(4,3); 


% loop for block
for jblock=block1:blockn

	% load file
	lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	load(lfile);


    % number of simulations in each block
	
    for jj=1:nsimul
        para=parasim(100*jj,:); 
alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a  = para(8);
rho_f  = para(9);
rho_i  = para(10);
sd_a1   = para(11);
sd_f1   = para(12);
sd_i1   = para(13);
sd_a2   = para(14);
sd_f2   = para(15);
sd_i2   = para(16);
piess  = para(17);
gy     = para(18);
lnA_0  = para(19);
yst    = para(20);
p11 = para(21);
p22 = para(22);

    
P_i = [p11 1-p11;
       1-p22 p22];
   

%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------

  A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
        -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];

    solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

    A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
        -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_f = inv(A_Af_mat) * [0;0;-kappa1;-kappa1;0;0];

    A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
        -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];

   
    y1a  = solvec_a(1);   y2a = solvec_a(2);
    y1f  = solvec_f(1);   y2f = solvec_f(2);
    y1i  = solvec_i(1);   y2i = solvec_i(2); 
    pi1a = solvec_a(3);  pi2a = solvec_a(4);
    pi1f = solvec_f(3);  pi2f = solvec_f(4); 
    pi1i = solvec_i(3);  pi2i = solvec_i(4);
    i1a  = solvec_a(5);   i2a = solvec_a(6);
    i1f  = solvec_f(5);   i2f = solvec_f(6);
    i1i  = solvec_i(5);   i2i = solvec_i(6); 

    var_a1 = (sd_a1^2)/(1-rho_a^2);   var_a2 = (sd_a2^2)/(1-rho_a^2);
    var_u1 = (sd_f1^2)/(1-rho_f^2);   var_u2 = (sd_f2^2)/(1-rho_f^2);
    var_e1 = (sd_i1^2)/(1-rho_i^2);   var_e2 = (sd_i2^2)/(1-rho_i^2); 


    Vout(1,:,(jblock-block1)*nsimul+jj)=[y1a^2*var_a1 y1f^2*var_u1 y1i^2*var_e1]./(sum(y1a^2*var_a1+y1f^2*var_u1+y1i^2*var_e1));
    Vout(2,:,(jblock-block1)*nsimul+jj)=[y1a^2*var_a2 y1f^2*var_u2 y1i^2*var_e2]./(sum(y1a^2*var_a2+y1f^2*var_u2+y1i^2*var_e2)); 
    Vout(3,:,(jblock-block1)*nsimul+jj)=[y2a^2*var_a1 y2f^2*var_u1 y2i^2*var_e1]./(sum(y2a^2*var_a1+y2f^2*var_u1+y2i^2*var_e1));
    Vout(4,:,(jblock-block1)*nsimul+jj)=[y2a^2*var_a2 y2f^2*var_u2 y2i^2*var_e2]./(sum(y2a^2*var_a2+y2f^2*var_u2+y2i^2*var_e2)); 
  
    Vinf(1,:,(jblock-block1)*nsimul+jj)=[pi1a^2*var_a1 pi1f^2*var_u1 pi1i^2*var_e1]./(sum(pi1a^2*var_a1+pi1f^2*var_u1+pi1i^2*var_e1));
    Vinf(2,:,(jblock-block1)*nsimul+jj)=[pi1a^2*var_a2 pi1f^2*var_u2 pi1i^2*var_e2]./(sum(pi1a^2*var_a2+pi1f^2*var_u2+pi1i^2*var_e2)); 
    Vinf(3,:,(jblock-block1)*nsimul+jj)=[pi2a^2*var_a1 pi2f^2*var_u1 pi2i^2*var_e1]./(sum(pi2a^2*var_a1+pi2f^2*var_u1+pi2i^2*var_e1));
    Vinf(4,:,(jblock-block1)*nsimul+jj)=[pi2a^2*var_a2 pi2f^2*var_u2 pi2i^2*var_e2]./(sum(pi2a^2*var_a2+pi2f^2*var_u2+pi2i^2*var_e2)); 
    
    Vint(1,:,(jblock-block1)*nsimul+jj)=[i1a^2*var_a1 i1f^2*var_u1 i1i^2*var_e1]./(sum(i1a^2*var_a1+i1f^2*var_u1+i1i^2*var_e1));
    Vint(2,:,(jblock-block1)*nsimul+jj)=[i1a^2*var_a2 i1f^2*var_u2 i1i^2*var_e2]./(sum(i1a^2*var_a2+i1f^2*var_u2+i1i^2*var_e2)); 
    Vint(3,:,(jblock-block1)*nsimul+jj)=[i2a^2*var_a1 i2f^2*var_u1 i2i^2*var_e1]./(sum(i2a^2*var_a1+i2f^2*var_u1+i2i^2*var_e1));
    Vint(4,:,(jblock-block1)*nsimul+jj)=[i2a^2*var_a2 i2f^2*var_u2 i2i^2*var_e2]./(sum(i2a^2*var_a2+i2f^2*var_u2+i2i^2*var_e2)); 
    

    end  
end
    doutmean(1,1)=mean(Vout(1,1,:)); doutmean(1,2)=mean(Vout(1,2,:)); doutmean(1,3)=mean(Vout(1,3,:));
    doutmean(2,1)=mean(Vout(2,1,:)); doutmean(2,2)=mean(Vout(2,2,:)); doutmean(2,3)=mean(Vout(2,3,:));
    doutmean(3,1)=mean(Vout(3,1,:)); doutmean(3,2)=mean(Vout(3,2,:)); doutmean(3,3)=mean(Vout(3,3,:));
    doutmean(4,1)=mean(Vout(4,1,:)); doutmean(4,2)=mean(Vout(4,2,:)); doutmean(4,3)=mean(Vout(4,3,:));
    
    drawout(:,1,1)=hpdint(squeeze(Vout(1,1,:)),hpdprob); drawout(:,1,2)=hpdint(squeeze(Vout(1,2,:)),hpdprob); drawout(:,1,3)=hpdint(squeeze(Vout(1,3,:)),hpdprob); 
    drawout(:,2,1)=hpdint(squeeze(Vout(2,1,:)),hpdprob); drawout(:,2,2)=hpdint(squeeze(Vout(2,2,:)),hpdprob); drawout(:,2,3)=hpdint(squeeze(Vout(2,3,:)),hpdprob); 
    drawout(:,3,1)=hpdint(squeeze(Vout(3,1,:)),hpdprob); drawout(:,3,2)=hpdint(squeeze(Vout(3,2,:)),hpdprob); drawout(:,3,3)=hpdint(squeeze(Vout(3,3,:)),hpdprob); 
    drawout(:,4,1)=hpdint(squeeze(Vout(4,1,:)),hpdprob); drawout(:,4,2)=hpdint(squeeze(Vout(4,2,:)),hpdprob); drawout(:,4,3)=hpdint(squeeze(Vout(4,3,:)),hpdprob); 
    
    dinfmean(1,1)=mean(Vinf(1,1,:)); dinfmean(1,2)=mean(Vinf(1,2,:)); dinfmean(1,3)=mean(Vinf(1,3,:));
    dinfmean(2,1)=mean(Vinf(2,1,:)); dinfmean(2,2)=mean(Vinf(2,2,:)); dinfmean(2,3)=mean(Vinf(2,3,:));
    dinfmean(3,1)=mean(Vinf(3,1,:)); dinfmean(3,2)=mean(Vinf(3,2,:)); dinfmean(3,3)=mean(Vinf(3,3,:));
    dinfmean(4,1)=mean(Vinf(4,1,:)); dinfmean(4,2)=mean(Vinf(4,2,:)); dinfmean(4,3)=mean(Vinf(4,3,:));
    
    drawinf(:,1,1)=hpdint(squeeze(Vinf(1,1,:)),hpdprob); drawinf(:,1,2)=hpdint(squeeze(Vinf(1,2,:)),hpdprob); drawinf(:,1,3)=hpdint(squeeze(Vinf(1,3,:)),hpdprob); 
    drawinf(:,2,1)=hpdint(squeeze(Vinf(2,1,:)),hpdprob); drawinf(:,2,2)=hpdint(squeeze(Vinf(2,2,:)),hpdprob); drawinf(:,2,3)=hpdint(squeeze(Vinf(2,3,:)),hpdprob); 
    drawinf(:,3,1)=hpdint(squeeze(Vinf(3,1,:)),hpdprob); drawinf(:,3,2)=hpdint(squeeze(Vinf(3,2,:)),hpdprob); drawinf(:,3,3)=hpdint(squeeze(Vinf(3,3,:)),hpdprob); 
    drawinf(:,4,1)=hpdint(squeeze(Vinf(4,1,:)),hpdprob); drawinf(:,4,2)=hpdint(squeeze(Vinf(4,2,:)),hpdprob); drawinf(:,4,3)=hpdint(squeeze(Vinf(4,3,:)),hpdprob); 
     
    dintmean(1,1)=mean(Vint(1,1,:)); dintmean(1,2)=mean(Vint(1,2,:)); dintmean(1,3)=mean(Vint(1,3,:));
    dintmean(2,1)=mean(Vint(2,1,:)); dintmean(2,2)=mean(Vint(2,2,:)); dintmean(2,3)=mean(Vint(2,3,:));
    dintmean(3,1)=mean(Vint(3,1,:)); dintmean(3,2)=mean(Vint(3,2,:)); dintmean(3,3)=mean(Vint(3,3,:));
    dintmean(4,1)=mean(Vint(4,1,:)); dintmean(4,2)=mean(Vint(4,2,:)); dintmean(4,3)=mean(Vint(4,3,:));
    
    
    drawint(:,1,1)=hpdint(squeeze(Vint(1,1,:)),hpdprob); drawint(:,1,2)=hpdint(squeeze(Vint(1,2,:)),hpdprob); drawint(:,1,3)=hpdint(squeeze(Vint(1,3,:)),hpdprob); 
    drawint(:,2,1)=hpdint(squeeze(Vint(2,1,:)),hpdprob); drawint(:,2,2)=hpdint(squeeze(Vint(2,2,:)),hpdprob); drawint(:,2,3)=hpdint(squeeze(Vint(2,3,:)),hpdprob); 
    drawint(:,3,1)=hpdint(squeeze(Vint(3,1,:)),hpdprob); drawint(:,3,2)=hpdint(squeeze(Vint(3,2,:)),hpdprob); drawint(:,3,3)=hpdint(squeeze(Vint(3,3,:)),hpdprob); 
    drawint(:,4,1)=hpdint(squeeze(Vint(4,1,:)),hpdprob); drawint(:,4,2)=hpdint(squeeze(Vint(4,2,:)),hpdprob); drawint(:,4,3)=hpdint(squeeze(Vint(4,3,:)),hpdprob); 
  xx=1:1:4;  
subplot(3,3,1); plot(xx,[doutmean(:,1) drawout(1,:,1)' drawout(2,:,1)']);     
subplot(3,3,4); plot(xx,[doutmean(:,2) drawout(1,:,2)' drawout(2,:,2)']);
subplot(3,3,7); plot(xx,[doutmean(:,3) drawout(1,:,3)' drawout(2,:,3)']); 

subplot(3,3,2); plot(xx,[dinfmean(:,1) drawinf(1,:,1)' drawinf(2,:,1)']);     
subplot(3,3,5); plot(xx,[dinfmean(:,2) drawinf(1,:,2)' drawinf(2,:,2)']);
subplot(3,3,8); plot(xx,[dinfmean(:,3) drawinf(1,:,3)' drawinf(2,:,3)']); 

subplot(3,3,3); plot(xx,[dintmean(:,1) drawint(1,:,1)' drawint(2,:,1)']);     
subplot(3,3,6); plot(xx,[dintmean(:,2) drawint(1,:,2)' drawint(2,:,2)']);
subplot(3,3,9); plot(xx,[dintmean(:,3) drawint(1,:,3)' drawint(2,:,3)']); 


